'data.frame': 26 obs. of 3 variables:
$ seedlings: int 7 10 0 5 54 0 30 14 1 2 ...
$ dungs : int 5 0 330 960 115 2330 475 10 2000 130 ...
$ property : Factor w/ 2 levels "Campesino","Empresa Forestal": 2 2 2 2 2 2 2 2 2 1 ...
Definición
Estructura
2.1. Función de vínculo
2.2. Familia de distribución de errores
Supuestos modelos lineales:
Los errores se distribuyen normalmente
La varianza (residual) es constante
La variable respuesta se relaciona linealmente con la(s) variable(s) independiente(s), si alguna de estas variables es continua
Uno o varios supuestos no se cumplen.
¿Qué hacer?
GLMs: extensión modelos lineales especificando estructuras no normales de errores y varianzas no constantes.
Modelos lineales caso particular de GLMs con distribución de errores normal y varianza constante.
Tipos de variables que por definición violan supuestos modelos lineales:
conteo de casos (p.ej. abundancia de especies)
conteo de casos como proporciones (p.ej. porcentaje plántulas muertas en un experimento)
respuesta binaria (p.ej. usado o no usado, presente o ausente)
3 componentes principales:
1. El predictor lineal
2. La función de vínculo
3. La estructura de los errores
Estructura modelos lineales:
yi = β0 + β1x1 + … + β1xn + Ɛi
Ɛi ~N(0,σ2)
Estructura GLMs:
f(yi) = β0 + β1x1 + … + β1xn + Ɛi
Var(y) ~ f(µ)
Transformación de la variable y buscando linealizar relación y con x y/o transformando en continuas variables que no lo son.
Por tanto, cambia la naturaleza de la variable respuesta.
Determinadas funciones de vínculo deben implementarse con determinados tipos de variables respuesta
Poisson: varianza aumenta linealmente con la media. Para respuestas tipo conteo.
Binomial negativa: extensión de Poisson. Para variables conteo con sobredispersión.
Binomial: para proporciones y datos presencia/ausencia.
Gamma: para datos continuos con coeficiente de variación constante.
Función glm()
Argumento family:
Estudio del efecto de la ganadería en la regeneración de araucaria. Muestreo de número de plantas de araucaria en dos tipos de ganadería: pequeños propietarios, con uso más intensivo de la tierra, y grandes empresas forestales, que mantienen algunos remanente de bosque nativo. Actividad ganado medida por número de bostas.
Variables:
Número de plántulas de araucaria
Número de bostas (proxy actividad del ganado)
Dos tipos de ganadería: pequeños propietarios y grandes empresas
Hipótesis:
H0: β = 0 (pendiente entre plántulas y número de bostas es 0)
H0: µ1 = µ2 (número plántulas en pequeños propietarios y empresas sin diferencias)
H0: β1 = β2 (relación entre plántulas y bostas es igual en pequeños propietarios y empresas)
glm_arau1 <- glm(seedlings ~ dungs*property, ara, family = "poisson")
anova(glm_arau1, test = "Chi")Analysis of Deviance Table
Model: poisson, link: log
Response: seedlings
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(>Chi)
NULL 25 347.56
dungs 1 30.967 24 316.59 2.624e-08 ***
property 1 80.598 23 236.00 < 2.2e-16 ***
dungs:property 1 37.974 22 198.02 7.170e-10 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Devianza explicada por el modelo
Medida de cantidad de variabilidad que explica el modelo (similar R2 en modelos lineales)
D2 = 1 - Dmodelo/Dnulo
Call:
glm(formula = seedlings ~ dungs * property, family = "poisson",
data = ara)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.497423 0.135408 18.444 < 2e-16 ***
dungs -0.010211 0.001789 -5.706 1.15e-08 ***
propertyEmpresa Forestal 0.576261 0.170467 3.380 0.000724 ***
dungs:propertyEmpresa Forestal 0.009041 0.001803 5.015 5.31e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for poisson family taken to be 1)
Null deviance: 347.56 on 25 degrees of freedom
Residual deviance: 198.02 on 22 degrees of freedom
AIC: 282.3
Number of Fisher Scoring iterations: 5
Fórmula: log(y) = 2.497 + 0.576 Empresa + (-0.010 + 0.009 Empresa) Bostas
Call:
glm.nb(formula = seedlings ~ dungs * property, data = ara, init.theta = 1.13715532,
link = log)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) 2.152593 0.362295 5.942 2.82e-09 ***
dungs -0.005663 0.002468 -2.294 0.0218 *
propertyEmpresa Forestal 1.113371 0.557993 1.995 0.0460 *
dungs:propertyEmpresa Forestal 0.003945 0.002528 1.561 0.1186
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for Negative Binomial(1.1372) family taken to be 1)
Null deviance: 49.506 on 25 degrees of freedom
Residual deviance: 28.040 on 22 degrees of freedom
AIC: 155.22
Number of Fisher Scoring iterations: 1
Theta: 1.137
Std. Err.: 0.384
2 x log-likelihood: -145.218
2 variantes de estos modelos:
Variable con valores binarios (1/0): modelos de regresión logística
Variables con valores porcentuales (nº éxitos / nº de ensayos)
Ejemplos:
Función de vínculo:
logit: log(p/(1-p))
logit(y) = ln(p/(1-p)) = β0 + β1x1 + … + β1xn
Curva sigmoidea entre 0-1
Estrategia reproductora de un alga invasora con dos especies nativas.
¿Probabilidad de encontrar estructuras reproductoras es diferente entre la especie invasora y las nativas teniendo en cuenta la talla del individuo?
Hipótesis: a misma talla del individuo, la especie invasora tendrá mayor probabilidad de reproducirse que las nativas
Variable respuesta (Estado):
- 0: sin estructuras reproductoras
- 1: con estructuras reproductoras
Variables predictoras:
- talla del individuo (cm)
- especie (Sp)
Analysis of Deviance Table
Model: binomial, link: logit
Response: Estado
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(>Chi)
NULL 358 476.38
Long 1 103.010 357 373.38 < 2.2e-16 ***
Sp 2 31.771 355 341.60 1.262e-07 ***
Long:Sp 2 0.459 353 341.15 0.7949
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Interacción no significativa
Analysis of Deviance Table
Model: binomial, link: logit
Response: Estado
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(>Chi)
NULL 358 476.38
Long 1 103.010 357 373.38 < 2.2e-16 ***
Sp 2 31.771 355 341.60 1.262e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
[1] 0.2829254
Call:
glm(formula = Estado ~ Long + Sp, family = "binomial", data = algas)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -2.80703 0.42701 -6.574 4.91e-11 ***
Long 0.26505 0.03174 8.351 < 2e-16 ***
SpTom -0.57476 0.34528 -1.665 0.096 .
SpVer -1.84877 0.36933 -5.006 5.57e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 476.38 on 358 degrees of freedom
Residual deviance: 341.60 on 355 degrees of freedom
AIC: 349.6
Number of Fisher Scoring iterations: 5
x_seq <- seq(min(algas$Long), max(algas$Long), length.out = 100)
newdata <- expand.grid(Long = x_seq, Sp = levels(algas$Sp))
# Obtener predicciones en escala de respuesta (probabilidades)
pred <- predict(glm_algas2, newdata, type = "link", se.fit = TRUE)
# Calcular intervalos de confianza
crit <- qnorm(0.975) # 95%
newdata$fit_link <- pred$fit
newdata$lwr_link <- pred$fit - crit * pred$se.fit
newdata$upr_link <- pred$fit + crit * pred$se.fit
# Convertir a escala de probabilidad (logística)
inv_logit <- function(x) exp(x) / (1 + exp(x))
newdata$fit <- inv_logit(newdata$fit_link)
newdata$lwr <- inv_logit(newdata$lwr_link)
newdata$upr <- inv_logit(newdata$upr_link)
# Promediar sobre el factor
pred_marginal <- aggregate(cbind(fit, lwr, upr) ~ Long, data = newdata, mean)
# Graficar relación promedio con IC
library(ggplot2)
ggplot(pred_marginal, aes(x = Long, y = fit)) +
geom_ribbon(aes(ymin = lwr, ymax = upr), fill = "#0072B2", alpha = 0.2) +
geom_line(color = "#0072B2", linewidth = 1.3) +
labs(
x = "Talla individuo (cm)",
y = "Predicción probabilidad estructuras reproductoras",
) +
theme_minimal(base_size = 14)newdat <- expand.grid(
Long = seq(min(algas$Long), max(algas$Long), length = 100),
Sp = levels(algas$Sp)
)
# Predicciones de probabilidad
pred <- predict(glm_algas2, newdat, type = "link", se.fit = TRUE)
# Convertir de escala logit a probabilidad
newdat$fit <- plogis(pred$fit)
newdat$lwr <- plogis(pred$fit - 1.96 * pred$se.fit)
newdat$upr <- plogis(pred$fit + 1.96 * pred$se.fit)
ggplot(newdat, aes(x = Long, y = fit, color = Sp, fill = Sp)) +
geom_line(linewidth = 1.2) +
geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.15, color = NA) +
labs(
x = "Talla individuo (cm)",
y = "Predicción probabilidad estructuras reproductoras",
) +
theme_minimal(base_size = 13)Comparaciones entre especies
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: glm(formula = Estado ~ Long + Sp, family = "binomial", data = algas)
Linear Hypotheses:
Estimate Std. Error z value Pr(>|z|)
Tom - Fra == 0 -0.5748 0.3453 -1.665 0.217824
Ver - Fra == 0 -1.8488 0.3693 -5.006 < 1e-04 ***
Ver - Tom == 0 -1.2740 0.3126 -4.076 0.000129 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)